﻿#include<iostream>
#include<cmath>
#include<iomanip>
#include<fstream>
#include<sstream> 
#include<cstdlib>
#include <time.h>
#include <windows.h>
#include <sys/timeb.h>
#include <time.h>

using namespace std;

/*定义常量
 [dim] 维度？2维or3维
 [dir] 离散速度个数
 [nx] x方向上的网格数
 [ny] y方向上的网格数
 [growth] 布置固相随机生长核概率
 [porosity] 目标孔隙率
 [probability] 每个方向上的生长概率
 [e] 离散速度方向
 [sum] 统计固体点的个数
 [ran1] 每个格点生成核的概率
 [ran2] 每个离散速度方向上生长概率
*/
const int dim = 2;
const int dir = 9;
const int nx=799;  
const int ny=293;
double porosity = 0.7;
double growth = 0.0025;
double Pc = 0.005;
double probability[dir - 1] = { 
  5*Pc     ,Pc     ,5*Pc     ,Pc    ,
  Pc / 4.,Pc / 4.,Pc / 4.,Pc / 4.
};
int e[dir][dim] = { 
  {0,0},
  {1,0},{0,1},{-1,0},{0,-1},
  {1,1},{-1,1},{-1,-1},{1,-1} 
};

int i,j,k,n;
double m,p,sum;
double ran1[nx+1][ny+1], ran2[nx+1][ny+1][9];
int s[nx+1][ny+1],szh[nx+1][ny+1];


void init();
void output();
double porous();
void character();

void init()
{
	sum=0;
	for(j=0;j<=ny;j++)
		for(i=0;i<=nx;i++)
		{
			s[i][j]=0;
			szh[i][j]=0;
		}
}

int main()
{
	void init();

	srand(time(0));
	for(j=0;j<=ny;j++)
		for(i=0;i<=nx;i++)
		{
			ran1[i][j]=(rand()%1000)/1000.0;
      if (ran1[i][j] <= growth)
      {
        s[i][j] = 1;
        szh[i][j] = 1;
        sum = sum + 1;
      }
		}

	sum=sum/(nx+1)/(ny+1);
	sum=1-sum;

	cout<<"生长核数"<<sum<<endl;

	ostringstream name1;
	name1<<"生长核.dat";
	ofstream out(name1.str().c_str());
	out<<"variables=X,Y,S"<<endl;
	out<<"zone I="<<nx+1<<"  "<<" J="<<ny+1<<endl;

	for(j=0;j<=ny;j++)
		for(i=0;i<=nx;i++)
		{
			out<<i<<"     "<<j<<"     "<<s[i][j]<<endl;
		}

  /*
  晶体生长
  */
  while (sum >= porosity)
  {
    srand(time(0));
    for (j = 0; j <= ny; j++)
      for (i = 0; i <= nx; i++)
        if (szh[i][j] == 1)
        {
          for (k = 1; k <= 8; k++)
          {
            ran2[i][j][k] = (rand() % 10000) / 10000.0;
          }
        }

    for (j = 0; j <= ny; j++)
      for (i = 0; i <= nx; i++)
      {
        if (szh[i][j] == 1)
        {
          for (k = 1; k <= 8; k++)
          {
            p = ran2[i][j][k];
            if (p < probability[k - 1])
              s[i + e[k][0]][j + e[k][1]] = 1;
          }
        }
      }

    porous();

    for (j = 0; j <= ny; j++)
      for (i = 0; i <= nx; i++)
      {
        szh[i][j] = s[i][j];
      }
  }

  /*
  过滤格点的信息
  */
	for(m=1;m<=2;m++)
	{
	for(j=1;j<ny;j++)
		for(i=1;i<nx;i++)
		{
			int sh[4];
			for (k=1;k<4;k++)
			{
				sh[k]=s[i+k][j]+s[i][j+k]+s[i-k][j]+s[i][j-k]+s[i+k][j+k]+s[i-k][j+k]+s[i-k][j-k]+s[i+k][j-k];
			}
			//int sh=s[i+1][j]+s[i][j+1]+s[i-1][j]+s[i][j-1]+s[i+1][j+1]+s[i-1][j+1]+s[i-1][j-1]+s[i+1][j-1];
			int sh1=s[i+1][j]+s[i][j+1]+s[i-1][j]+s[i][j-1];
			if(s[i][j]==1)
			{
				if(sh[1]<=1)
				{
					s[i][j]=0;
				}
			}
			else
			{
				if(sh[1]>=7)
				{
					s[i][j]=1;
				}
				if(sh1==4)
				{
					s[i][j]=1;
				}
				for(k=1;k<4;k++)
				{
					if(sh[k]==8)
					{s[i][j]=1;}
				}
			}
		}
	}
	output();
	character();
	system("pause");
	return 0;
}

void output()
{
	ostringstream name2;
	name2<<"2D.dat";
	ofstream out(name2.str().c_str());
	out<<"variables=X,Y,S"<<endl;
	out<<"zone I="<<nx+1<<"  "<<" J="<<ny+1<<endl;

	for(j=0;j<=ny;j++)
		for(i=0;i<=nx;i++)
		{
			out<<i<<"     "<<j<<"     "<<s[i][j]<<endl;
		}

  ostringstream name;
  name << "matlab.dat";
  ofstream out1(name.str().c_str());
  for (j = 0; j <= ny; j++)
    for (i = 0; i <= nx; i++)
    {
      out1 << i << "     " << j << "     " << s[i][j] << endl;
    }
}

double porous()
{
		sum=0;
		for(j=0;j<=ny;j++)
			for(i=0;i<=nx;i++)
			{
				if(s[i][j]==1)
				{
					sum=sum+1;
				}				
			}

		sum=sum/(nx+1)/(ny+1);
		sum=1-sum;
		return sum;
}

void character()
{
	sum=porous();
	ostringstream name3;
	name3<<"porosity.dat";
	ofstream out(name3.str().c_str());
	out<<"porosity="<<sum<<endl;
}
